function [ A, Z ] = assemblyMatrix( domain2D )
%ASSEMBLYMATRIX Summary of this function goes here
%   Detailed explanation goes here

% nV = size(domain2D.vx, 1);
nE = size(domain2D.e1, 1);

A = diag(ones(2*nE, 1)) * 0.5;
Z = ones(2*nE, 2*nE);

for i = 1 : nE
    
   pi0 = domain2D.e1(i);
   pi1 = domain2D.e2(i);
   
   p0 = [domain2D.vx(pi0), domain2D.vy(pi0)];
   p1 = [domain2D.vx(pi1), domain2D.vy(pi1)];
   
   p = (p0 + p1) / 2.0;
    
%    p = [domain2D.vx(i), domain2D.vy(i)];
   
   for j = 1 : nE
      qi0 = domain2D.e1(j);
      qi1 = domain2D.e2(j);
      
      q0 = [domain2D.vx(qi0), domain2D.vy(qi0)];
      q1 = [domain2D.vx(qi1), domain2D.vy(qi1)];
      
      qn = [domain2D.nx(j), domain2D.ny(j)];
      
      v0 = evaluateGradientQuadraturePoint(p, q0, qn);
      v1 = evaluateGradientQuadraturePoint(p, q1, qn);
      
      ql = domain2D.l(j);
      
      A(2*i-1:2*i, 2*j-1:2*j) = A(2*i-1:2*i, 2*j-1:2*j) + (v0 + v1) * ql / 2;
      
   end
   
%    for j = 1 : n
%        qi0 = j;
%        qi1 = mod(j, n) + 1;
%        
%        q0 = [domain2D.x(qi0), domain2D.y(qi0)];
%        q1 = [domain2D.x(qi1), domain2D.y(qi1)];
%        
%        qn = [domain2D.nx(j), domain2D.ny(j)];
%        
%        v0 = evaluateGradientQuadraturePoint(p, q0, qn);
%        v1 = evaluateGradientQuadraturePoint(p, q1, qn);
%        
%        ql = domain2D.l(j);
%        
%        A(2*i-1:2*i, 2*j-1:2*j) = A(2*i-1:2*i, 2*j-1:2*j) + (v0 + v1) * ql / 2;
%        
%    end
   
   
   
end

A = -1 * A;

end

